#include<bits/stdc++.h>

using namespace std;
const int MAX=2e3+3;
const double min_rock_d=1.5;
const double min_crater_d=10;
const int length=2e3;
const double resolution=0.5;
const double smoothness=25;
const double h1_d[3][2]= {0.23,0.17,0.11,0.25,0.19,0.13}; //h1_d[i][0] for lower_limit, h1_d[i][1] for upper_limit
const double h2_d[3][2]= {0.022,0.016,0.008,0.06,0.045,0.03}; // ditto


double dem[MAX][MAX];
bool vis[MAX][MAX];
struct node
{
    int xs,xt,ys,yt;
    double s;
    node(int _xs,int _xt,int _ys,int _yt,double _s):xs(_xs),xt(_xt),ys(_ys),yt(_yt),s(_s){}
};

inline double sq(double x)
{
    return x*x;
}

inline double rock(double d,double h,double x,double y)
{
    double res=h-4.0*h*(sq(x)+sq(y))/sq(d);
    if(res<0)return 0.0;
    return res;
}

inline double crater(double d,double h1,double h2,double x,double y)
{
    x=sqrt(sq(x)+sq(y));
    double res;
    double rd=rand()%2+1;
    if(x+rd<=d/2)
        res=4.0*h1/sq(d)*sq(x)+h2-h1;
    else
        res=-4.0*h2/sq(d)*sq(x)+4.0*h2/d*(x);
    if(x+rd>d/2&&res<0)
        return 0.0;
    return res;
}

inline int c_d_n(double d)
{
    if(d<=40.0)
        return int(pow(10.0,-2.0*log10(d)-1.0)*length*length/4);
    else
        return int(pow(10.0,-3.0*log10(d)+0.62)*length*length/4);
}
inline double c_n_d(int n)
{
    double nn=double(n)/sq(double(length/2));
    if(n<65)
        return pow(10.0,(log10(nn)-0.62)/-3.0);
    else
        return pow(10.0,(log10(nn)+1.0)/-2.0);
}
inline int r_d_n(double d)
{
    return int(pow(10.0,-2.589*log10(d)-3.3955)*length*length/4);
}
inline double r_n_d(int n)
{
    double nn=double(n)/sq(double(length/2));
    return pow(10.0,(log10(nn)+3.3955)/-2.589);
}

void fill_rock(double d,double h,int xs,int xt,int ys,int yt,int xm,int ym)
{
    if(xs<0)xs=0;
    if(xt>=length)xt=length-1;
    if(ys<0)ys=0;
    if(yt>=length)yt=length-1;
    for(int i=xs; i<=xt; i++)
        for(int j=ys; j<=yt; j++)
        {
            double xx=double(i-xm)*resolution,yy=double(j-ym)*resolution;
            dem[i][j]+=rock(d,h,xx,yy);
        }
}
void fill_crater(double d,double h1,double h2,int xs,int xt,int ys,int yt,int xm,int ym)
{
    /*
        d:   坑的直径
        h1： 坑内高度
        h2： 边缘高度
        xs： x-r
        xt： x+r
        ys： y-r
        yt： y+r
        xm： x
        ym： y
    */
    if(xs<0)xs=0;
    if(xt>=length)xt=length-1;
    if(ys<0)ys=0;
    if(yt>=length)yt=length-1;
    for(int i=xs; i<=xt; i++)
        for(int j=ys; j<=yt; j++)
        {
            double xx=double(i-xm)*resolution,yy=double(j-ym)*resolution;
            dem[i][j]+=crater(d,h1,h2,xx,yy);
            vis[i][j]=1;
        }
}

void terrain_genaration(void)
{
    srand(time(0));
    queue<node>Q;
    Q.push(node(0,length-1,0,length-1,smoothness));
    dem[0][0]=dem[0][length-1]=dem[length-1][0]=dem[length-1][length-1]=10.0;
    while(!Q.empty())
    {
        node cur=Q.front();
        Q.pop();
        if(cur.xt-cur.xs<2)continue;
        int xs=cur.xs,xt=cur.xt,ys=cur.ys,yt=cur.yt;
        int mx=(xs+xt)>>1,my=(ys+yt)>>1;
        double rd=double(rand()%int(max(cur.s*2000.0,1.0)))/1000.0-cur.s;
        double avg=(dem[xs][ys]+dem[xs][yt]+dem[xt][ys]+dem[xt][yt])/4.0;
        dem[mx][my]=avg+rd;
        rd=double(rand()%int(max(cur.s*2000.0,1.0)))/1000.0-cur.s;
        dem[xs][my]=avg+rd;
        rd=double(rand()%int(max(cur.s*2000.0,1.0)))/1000.0-cur.s;
        dem[xt][my]=avg+rd;
        rd=double(rand()%int(max(cur.s*2000.0,1.0)))/1000.0-cur.s;
        dem[mx][ys]=avg+rd;
        rd=double(rand()%int(max(cur.s*2000.0,1.0)))/1000.0-cur.s;
        dem[mx][yt]=avg+rd;
        double nxts=cur.s/2.0;
        Q.push(node(xs,mx,ys,my,nxts));
        Q.push(node(mx,xt,ys,my,nxts));
        Q.push(node(xs,mx,my,yt,nxts));
        Q.push(node(mx,xt,my,yt,nxts));
    }
}

void Creater(int noc,int a,int b)
{
    double minn=1e9+7;
    for(int i=0;i<a;i++)
    {
        int num=rand()%noc+b;
        double d=c_n_d(num);
        minn=min(minn,d);
        int t=rand()%3;
        double h1=h1_d[t][0]+(h1_d[t][1]-h1_d[t][0])*double(rand()%10+1)/10.0;
        h1*=d;
        double h2=h2_d[t][0]+(h2_d[t][1]-h2_d[t][0])*double(rand()%10+1)/10.0;
        h2*=d;
        int r=int(d/resolution);
        int x=rand()%length,y=rand()%length;
        while(vis[x][y])
        {
            x=rand()%length;
            y=rand()%length;
        }
        fill_crater(d,h1,h2,x-r,x+r,y-r,y+r,x,y);
    }
}

int main()
{
    srand(time(0));
    terrain_genaration();         // diamond-square algorithm
    int noc=c_d_n(min_crater_d);
    cout<<"noc="<<noc<<endl;
    Creater(noc,noc*0.3,noc*2.5);
    Creater(noc,noc*0.2,noc*0.5);
    Creater(noc,noc*0.15,2);
    int nor=r_d_n(min_rock_d);
    cout<<"nor="<<nor<<endl;
    for(int i=0; i<nor; i++)
    {
        int num=rand()%nor+1;
        double d=r_n_d(num);
        double h=double(rand()%50+50)/100.0*d;
        int r=int(d/resolution)/2;
        int x=rand()%length,y=rand()%length;
        fill_rock(d,h,x-r,x+r,y-r,y+r,x,y);
    }
    FILE *fpwrite=fopen("maptest.txt","w");
    for(int i=0;i<500;i++)
        for(int j=0;j<500;j++)
        {
            int si=i*4,ti=i*4+4,sj=j*4,tj=j*4+4;
            for(int k=si;k<ti;k++)
                for(int l=sj;l<tj;l++)
                    fprintf(fpwrite,"%.5f %.5f %.5f\n",double(k)*0.5,double(l)*0.5,dem[k][l]);
        }
}


